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ABSTRACT 

This  paper  introduces  methods  tailored  especially  for  problems 

'■■  x 
whose  solution  behaves  like  e   ,  where  \    is  complex.   The  shallow 

water  equations  with  topography  admit  such  solution. 

This  paper  complements  the  results  of  Pratt  and  others  on 

exponential-fitted  methods  and  those  of  Gautschi,  Neta,  van  der 

Houwen  and  others  on  trigonometrically-f itted  methods. 

1 .   Introduction 

In  this  paper  we  consider  linear  multistep  methods 


a,y       =  h   J   bnf(x   -   ,y  .,  , ),  k    1,  n    k-1       (1) 
z  =  [)  i.  =  (J 


for  integrating  the  initial  value  problem 

y'  (x)  =  f  (x,y  (x)  )  ,   y(xQ)  =  yQ  .  (2 

This  linear  multistep  method  is  characterized  by  the  polynomials 


0(O    =      V   a„cW  ,    cC)  =     b-,-k"   .  (3) 


The  main  assumption  of  this  paper  is  that  it  is  a  priori  known 
that  the  solution  is  approximately  of  the  form 


m      i>  •  t 
v(x)  -  cn  +      c.e  (4) 


where  A.  =  w.  +  ity . ,  and  the  frequencies  w.  are  in  a  given  interval 
[w£,wu] . 

The  special  case  where  >  .  =  jwQ  with  w„  given  was  considered 
first  by  Gautschi  [8] .   His  approach  was  the  following.   Let: 


i  (z)  =  p  (eZ)  -  zc(eZ)  (5) 


then  the  local  truncation  error  of  (1)  is  given  by  Lambert  [11] 


Tn+k  "  '(h  Ht>  y(tn>  •  <6» 


Insertina  (4)  in  (6)  vields 


Tn+k  .  ;(0)c0  -  _|   c.^ihX.le1  i       ,  A.  =  jw0  .  (7 


The  coefficients  b.  are  chosen  in  such  a  way  that 

4>  (ihjwQ)  =  0  ,  j  =  0,1,.  .  .  ,q  ,  (8) 

for  the  largest  value  of  q  possible.   q  is  then  called  the  trigo- 
nometric order  of  the  method.   Gautschi  has  chosen  a.  such  that  the 

i 

methods  are  of  Adams  and  Stormer  type.   However,  these  methods  are 
sensitive. to  changes  in  the  frequency  w~ .   Neta  and  Ford  [13] 
developed  Nystrom  and  generalized  Milne-Simpson  type  methods.   These 
methods  showed  less  sensitivity  to  perturbation  in  w„  but  require 
the  eigenvalues  of  the  Jacobian  to  be  purely  imaainary.   Neta  [14] 
has  developed  families  of  backward  differentiation  methods  that 
overcome  the  above-mentioned  restriction.   Salzer  [17]  has  developed 


predictor-corrector  methods  based  on  trigonometric  polynomials. 

See  also  Steifel  and  Bettis  [18]  and  Bettis  [3] .   Van  der  Houwen 

and  Sommeijer  [10J  have  developed  an  alternative  approach.   The 
conditions  (8)  were  replaced  by 

<M0)  =0  , 

(9) 

<J>(ihA  (j)  )  =  0  ,   j  =  1,2,.  .  .  ,q  , 

where  the  A     are  appropriately  chosen  points  in  the  interval 
[w  o ,w  ]  . 

An  advantage  of  this  so-called  minimax  approach  over  the  fittinq 
approach  is  the  increased  accuracy  in  cases  where  no  accurate  esti- 
mate of  w~  is  available  or  when  the  frequency  is  varying  in  time. 
The  ether  special  case  considered  in  the  literature  is  where 
>.  .  =  iv  •  .   Probably  the  first  article  on  the  subject  is  due  to  Brock 
and  Murray  [5].   They  discuss  the  use  of  exponential  sums  in  the 
integration  of  a  system  of  first  order  ordinary  differential  equations 
Dennis  [7]  also  suggested  special  methods  for  problems  whose  solution 
is  exponential.   He  suggested  a  transformation  of  variables.   More 
recently,  Carroll  [6J  has  developed  exponentially  fitted  one-step 
methods  for  the  scalar  Riccati  equation.   For  the  general  first 
order  system  of  equations,  Pratt  [16]  suggests  methods  based  on  the 
three  parameter  exponential  function 


I (x)  =  A  +  BeZX  .  (10) 


The  parameters  A,  B  are  given  in  terms  of  values  of  y  and  f. 
Several  possibilities  for  z  are  given  based  on  results  of  Brandon 
[2]  and  Babcock  et  al .  [1]  . 

Lyche  [12]  analyzes  multistep  methods  which  exactly  integrate 

w  x 
.  .       r  m   n  i    ,         .      ,     . 
the  set  tx  e    },  where  w   is  real  or  imaainary. 

n  2 

In  this  article  we  developed  various  methods  fitting  exponentials 
and  methods  obtained  via  the  minimax  approach. 

2 .   Construction  of  Methods 
2 . 1  Fitting  Methods 

In  this  subsection  we  discuss  various  fitting  methods.   To  this 
end,  we  separate  MihjX)  =  0,  j  =  1,2,.. .,q,  into  real  and  imaginary 
parts.   This  yields  the  following  equations  relating  the  coefficients 
a. ,br  , 


aej~(k  £)cos  jv(k_o)  _   y  b,er,j(k  l)   [ju  cos  jv(k-£) 

v=o  •'■  e=o  k 

-    jv   sin  jv  (k-£) ]    =    0    , 


V  k 

j  y  ( k  -  £ )  , ,       r  .  r     ,  j  \i  ( k  - 1 )    r  .  .        .      n       „  N 

a,eJ  sm   ]v(k-£)    -       /    b0eJ  [jy   sin  jv  (k-£ ) 


(11) 


L 
£=0  1=0 


+    j v  cos  jv  (k-£) ]  =  0  , 


where  /  =  w  +  i^ ,  y  =  -h^ ,  v  =  hw ,   j  =  l,2,...,q. 

For  explicit  methods,  b   =  0.   For  Adams  type  methods  a_  =  1,  a,  =  -1 

a.  =  0  for  i  =  2,...,k.   For  Nystrom  or  generalized  Milne-Simpson 

methods  a-  =  1.  a-  =  -1  and  other  a.  =0. 
0        2  l 


k  =  1  Implicit 


Adams 


(12) 


h     ve    4-  u  sinv  -  v  cosv 
(u  +v  )  sin  v 


,     ve   -  u  smv  -  v  cosv 

bl  =  2 2 

(V  +  v  )  sin  v 


—  COS  \' 

For  iL  =  0,  the  coefficients  become  b0  =  , =  b,  ,  which  aaree 

0      sin  v       1 

with  Gautschi  [8]  if  the  coefficients  are  expanded  in  Taylor  series 
with  respect  to  v. 

k  =  2  Explicit 
Adams 

(u   sin2v    -    v  cos2v)e"    +    v  cosv    -    u   sinv 


b,    = 


1  ,22,. 

(y    +\     )     sin    v 


>    2l 

(v  cosv    -    u   sinv)e         -    ve 

2  2       ? 

(u    +v    )     sin    v 


(13) 


Nvstrom 


,        _     (u   sin2v    -    v   cos2v)eM    +    ve         , 

1  2      2 

(jj    +v    )     sin    v 


/  .        ,     2p  -u 

,  (v  cos  v    -    u   sinv)e         +    ve 

b~    = 


(14) 


2  2       2. 

(u    +v    )     sin    v 


sm    v 


For    'j,    =    0 ,    the   coefficients   become    b,    =    2    — ,    b2    =    0   which 

aqree   with    Neta     [14]  . 


k  =  2  Implicit 

In  this  case,  one  obtains  a  one-parameter  family  of  (Adams, 
generalized  Milne-Simpson)  methods  of  trigonometric  order  1.   The 
free  parameter  can  be  used  to  increase  the  algebraic  order  of  the 
method  as  in  [13] . 

Backward  Differentiation 

e    cos  2v  +  a, e   cos  v  +  a?    -  b„e   (y  cos2v  -  v  sin2v)  =  0  , 


e  '    sin  2v  +  a,e   sin  v       -  b~e  "  (u  sin2v  +  v  cos2v)  =  0  ,    (15 


+  a,  +  a  =  0 


This  system  can  be  solved  by  MACSYMA  (Project  MAC ' s  SYmbolic  MAnipu- 
lation  system  written  in  LISP  and  used  for  performing  symbolic  as 
well  as  numerical  mathematical  manipulation  [4])  or  by  REDUCE  [9]. 
The  solution  is 


ve    -  u  sin2v  -v  cos2v 
a   =  — —  , 

-e  (v  cosv  +  y  sinv)  +  p  sm2v  +  v  cos/v 

(16) 


■  -e    sin  v  +  e   sin2v  -  sin  v 

h   =  — — ■ — ■ 

-e    ^  (v  cosv    +    \i   sinv)    +    ep  (y   sin2v    +    v  cos2v) 


For    \p   =    0 ,    the    coefficients    agree   with    those    given    in     [14J 


k  =  3  Explicit 

Again  here,  one  obtains  a  one-parameter  family  of  methods  of 
trigonometric  order  1.   In  order  to  get  methods  of  trigonometric 
order  2,  one  has  to  construct  a  3  step  implicit  method  of  Adams  or 
generalized  Milne-Simpson  type.   In  order  to  increase  the  trigo- 
nometric order  without  going  to  a  higher  step  number,  one  can  con- 
struct linear  multistep  methods  for  which  the  coefficients  a.  are 
also  functions  of  ty ,    w.   Some  examples  are  given  in  the  next 
subsection . 

2 . 2   Generalized  Fitting  Methods 

In  this  section,  we  construct  some  linear  mulstistep  methods 
of  the  form 


k  k 

I   a;.  (A)y     .  =  h   [ 
-0        n  £=0 


a.U)yn+1_£  =  h  I   V>.)fn+1_,  ,   k  »  1,  n  >  k-1  .    (17) 


Since  a,  are  functions  of  X    one  has  more  free  parameters  for  his 
disposal  which  can  be  used  to  obtain  higher  trigonometric  order 
methods  with  relatively  lower  step  number. 

k  =  2  Implicit 

In  this  case,  one  has  to  solve  the  following  linear  system  of 
five  equations  for  the  parameters  a,,  a^  ,    bQ ,  b^ ,  b2  to  obtain  a 
method  of  trigonometric  order  2. 


1  +a1  +a2  =    0    , 

2y  -,  y  ,       2y  .  n  ~    . 

e        cos2v  +  a,e      cos  v  +  a~   ~"Ae       (y  cos2v   -v   sin2v) 

-b,e     (y  cosv   -v   sinv)    -  ub„    =    0    , 

2u  n       ,  P  t_       2y  .         .     _  .     . 

e        sm2v  +  a,e      sinv  ~^ne       ^u   sin2v  +v  cos2v) 

-b,e  '  (y   sinv  +  v  cosv)    -vb~    =    0    ,  (18 

4y  2u  4u 

e     '  cos4v  +a,e     '  cos  2v+a„- b»e       (2y   cos4v  -  2v   sin4v) 

2y 
-b,e       (2y  cos2v   -2v   sin2v)-2ybp   =0    , 

4  u  2  u  4 ! 

e        sin4v  +  a,e        sin2v        ~^ne    "  ^u   sin^v  +2v  cos4v) 


2u 


-b-,e   (2-M  sin2v  +2v  cos2v)-2vb0  =  0 


The  system  was  solved  by  REDUCE  [9].   The  expressions  for  the 
coefficients  are  complicated  but  REDUCE  produces  an  output  in  the 
form  of  Fortran  statements  that  can  be  incorporated  in  a  computer 
program  for  numerical  experiments  with  such  a  method. 

2  .  3   Mmimax  Methods 

In  this  section  we  discuss  minimax  methods,  i.e.,  methods 
obtained  by  satisfying  conditions  (9).   These  conditions  can  be 
written  in  terms  of  a0,  bp  as  follows: 

A,         A/ 

I    a,e;  '-)u         cos(k-£)vl]J 
£  =  0  ** 

b.e(k-£)v(D){-vJj,sin(k-Ov(j)+u(j,cos(k-Jl)v(j)},         (19) 
i.  =  0    ' 

j  =  1 , 2 , . . . ,q  , 


Z„«(k"1,,lt3)«in<lc-t>v<J> 


1=  l 


1  =  ^ 


k  ( i ) 

b0e(k"e)w    (v(j)cos(k-l)v(j)+li(j)sin(K)v(])  },       (20 

j  =  1 , 2  ,  .  .  .  ,q  , 


where  y(^  =  he  ^  )  ,    v(j)  =  hw ( 3 >  . 

ih,\  J   are  the  zeros  of  the  function  £(ih>)  such  that  it  has 
a  small  maximum  norm  in  the  rectangle  w.  <  w    wu ,  typ         .    • 

To  obtain  the  best  approximation  in  this  case  is  certainly 
not  easy;  but  we  will  assume  that  one  can  write  ?(ih\  -*  )  = 
(My1    >w    )  as  a  product  of  2  one  variable  functions.   Thus 
ty    3      and  w  -1   can  be  taken  as  Chebyshev's  points  on  the  corres- 
ponding interval,  i.e. 


(j)    1  ,       >    1  ,     ,  ,     2i-l 
J   =  T  -   +  .     +  7T  .'   -yp  cos  -^ — 
2  l         u     2   u    t       2q 


(21 


W  J    =  tt(W„  +W     +  Tr(VJ        -W,  COS  — k T   , 

2   I    u     2   u    -I       2q 

j  =  1,2  ,  .  .  .  ,q  . 

For  this  choice  of  points,  one  can  evaluate  the  coefficients 
a„  ,  b0.  by  solving 

c  (0)  =  o  , 

(22) 
:  (f  (D)  ,w(j))  =  0  ,   j  =  1,2,.  .  .  ,q  . 

We  call  such  methods  product  minimax  (PM*"). 
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The  number  of  free  parameters  for  implicit  methods  2k+l  and 
the  number  of  equations  is  2q+l,  thus,  the  trigonometric  order 
q  is  equal  to  the  step  number  k. 

k  =  1  Implicit 

In  this  case 


4>  =    j^  I   +^u) 


w(1)  =  ^(w£  +wu)  (23) 


a0  =  1  =  -ax 


and  the  system  of  equations  can  be  solved  for  b  ,  b,  .   This 
yields  the  coefficients  given  by  (12)  where 

m  =  -Yi%  (1)/v  =  hw(1).  (24 

Thus,  the  product  minimax  method  would  suggest  using  the  center 
of  the  rectangle  [^  n  ,i>    ]x[w„,w  ]  as  XQ.   To  obtain  a  product 
minimax  method  of  trigonometric  order  2 ,  one  has  to  solve  a 
system  of  5  equations  similar  to  equation  (18)  with  the  unknowns 

b„,  b, ,  b? ,  a, ,  a_  .   The  difference  is  that  in  the  last  2  equa- 

(2 )  (2  ) 

tions  one  should  replace  2y  by  u     and  2v  bv  v    .   In  the 

second  and  third  equations  of  (18),  the  u ,  v,  should  be  replaced 

by  p    ,  v     respectively.   The  resulting  system  can  be  solved 

by  MACSYMA  [4]  or  REDUCE  [9]. 

In  the  next  section  we  implement  two  methods  of  trigonometric 

order  1  and  2  and  see  how  the  product  minimax  methods  compare 

with  fitting  methods. 
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3 .   Numerical  Example 

Both  systems  (18)  and  (19)- (21)  were  solved  by  REDUCE  which 
produced  a  FORTRAN  subroutine  for  the  evaluation  of  the  coeffi- 
cients.  This  subroutine  is  called  only  once  during  the  integration 

The  methods  were  compared  for  the  solution  of  the  initial 
value  problem 


z-yi^+i)2^0  /    °ltl4  ' 


z(0)  =  1 


(25) 


whose  exact  solution 


ij(l+i)t 
z(t)  =  e  ,  (26) 


thus 


A  =  j  ( 1  +  i )  ,   V  =  "  2  /     w  =  2  '  ( 2  7 


In  order  to  avoid  complex  arithmetic,  we  rewrite  the  differential 
equation  as  a  system  of  equations  for  the  real  and  imaginary 
part  of  z  =  u  +  iv. 


u  +  jiu   + v)  =0  , 


0  <  t  <  4 


v  -  |(u  -v)  =  0  ,  (23 


u(0)  =  1  , 
v  (0)  =0  . 
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The  system  is  solved  by  fitting  methods  of  trigonometric  order 
1  and  2  with  h  =  .01  and  various  values  of  \p   and  w.   in  Table 
1  we  list  the  Euclidean  norm  of  the  error  at  t  =  4.   It  is 
clear  that  the  method  is  not  sensitive  to  perturbations  in  the 
values  of  i>    and  w. 


M> 

Aw 

error 

first  order 

second  order 

0 

0 

.3678  (-12) 

.4433  (-12) 

0 

.1 

.4482  (-7) 

.3508  (-11) 

.1 

0 

.4346  (-7) 

.2544  (-11) 

.1 

.1 

.6342  (-7) 

.442K-11) 

0 

.2 

.9241  (-7) 

.8126(-11) 

.2 

0 

.8706 (-7) 

.5095 (-11) 

Table  1 

Using  the  product  minimax  methods  of  trigonometric  order 
1  and  2  with  h  =  .  01  and  various  squares  centered  at  ip  =  -  tt/2  , 
w  =  tt/2,  the  error  is  much  laraer  but  again  is  insensitive 
to  small  perturbations  in  the  length  of  the  sides  of  the 
squares.   In  Table  2  we  list  the  Euclidean  norm  of  the  error 
at  t  =  4. 
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length  of  side 

error 

first  order 

second  order 

.4 

.3678(-12) 

.1469 (-6) 

.8 

.3678(-12) 

.2348 (-5) 

1.2 

.3678  (-12) 

. 1186  (-4) 

1.6 

.3678(-12) 

.  3728  (-4) 

2 

.3678(-12) 

.9001  (-4) 

Table  2 

Note  that  the  perturbations  in  the  product  minimax  methods  are 
larger  than  those  allowed  in  the  fitting  methods.   It  is 
possible  that  the  larger  errors  in  the  product  minimax  methods 
are  due  to  the  assumption  that  ;  can  be  v.ritten  as  a  product 

of  2  one  variable  functions.   Also  note  that  for  the  first 

2 
order  method,  one  always  gets  a  good  result  since  PM   always 

uses  the  center  of  the  square. 
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APPENDIX 


Here  we  show  that  the  shallow  water  equations  witn 
topography  have  a  solution  of  the  form  eAx,  where  A  is  complex. 
This  system  ot  equations  consists  of  three  equations  with 
three  forecast  variables  u,  v  and  $ .   The  equations  are: 


3u  3u  3u  3ff 

at  +  u  a3T  +  v  37  -  fv  +  g  -  °  .  (A.i) 


3v  3v  ?v  7)^ 

tt+Udt+V^+fu+~=°>  (A. 2) 


3(^  9     r      /  vi  3     r 

Ft    +    3x"Iu(^B)]     +    37[v(9--b)]     =    °     '  (A. 3) 


where  $  =  gh  is  the  geopotential  height  (h  =  height  of  free 
surface),  ^  is  the  bottom  topography  (assumed  to  be  independent 
of  time),  u,  v  are  the  components  of  the  wind  velocity  in  the 
x,  y  direction,  respectively,  and  f  is  the  Coriolis  parameter. 
Linearizing  the  equations  by  letting 

u  =  U  +  u  '  ,   v  =  V  +  v '  ,   c=  $+<$>'     , 

where  U,  V  are  the  constant  mean  flow  and  v  is  independent  of 
time.   Assuming  that  U,  V  are  related  to  $  via  the  geostrophic 
relations 


„  _  _  1  3*      ..    1  3$ 

U~    f  37  '    V=f3x"  (A. 4 


one  obtains  the  linear  system  (after  dropping  the  primes) : 


17 


at  dx      dy        9x      ' 


(A. 5) 


3t      dx      dy         dy 


0  , 


(A. 6) 


i±  +  u  |i  +  v  f*  +  y^ii  +  i^xi  =  U 

dt       d  X       dy      aX  dy  dX 


•'*  +  v  :Cb 


oy 


(A. 7) 


where  y   =    $  -  <p 


B 


If  the  flow  is  assumed  to  be  along  the  topography  as  in 
[15] ,  then  the  right  hand  side  of  (A. 7)  is  zero.  In  such  a 
case,  one  can  write  the  solution  in  the  form 


where 


u  =  u  e 
o 


v  =  v  e 
o 


i  (  E,  X  +  r  y  -at 


i  ( £ x  +  n  v  -at) 


ei  (?x  +ny  -c- 
*o 


I    =  V 


-    IP  , 


(A. 8) 


(A.  9) 


n  =  v  -  iG  . 

In  order  for  (A. 8)  to  be  a  solution  for  (A.5)-(A.7),  one  must 
have 


A  =  -o  +  ;U  +  nv 


(A. 10 


satisfying 


iA[(iA)2    -    in(inY    +  |J]    -    f[-iXf    -    iCdnY    +  I1)  J 

oy  dy 

+  (iCy  +  §J)  (-inf  +  ?A)  =  0  .  (A. 11) 


The  real  and  imaginary  parts  of 


A   =  -a  +  yU  +  vV  , 


(A. 12) 


a  .  =  -pu  -  ev  , 


satisfy  the  following  system  of  equations  (after  dropping 
nonlinear  terms  in  A) 


iU  dX  +  n  3yj    r  in  3x    *  3y;  ' 


Ar[f2  +  y(C2  +n2)]  +  A 

(A. 13 


X  (5  |X+n  |X)  -  a.  [f2  +y(?2  +n2)]  =  0  . 
r  '  dx  dv      i      '  ^ 


In  general,  >  is  complex  and,  thus,  the  shallow  water  equations 
have  a  solution  in  the  class  of  problems  to  be  discussed  here. 
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